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ABSTRACT 

We discuss early results from the first large N-body/hydrodynamical simulation to resolve the formation of 
galaxies in a cold dark matter universe. The simulation follows the formation of galaxies by gas cooling within 
dark halos of mass a few times IO'^Mq and above, in a flat universe with a positive cosmological constant. 
Over 2200 galaxies form in our simulated volume of (lOOMpc)^. Assigning luminosities to the model galaxies 
using a spectral population synthesis model results in a K-band luminosity function in excellent agreement with 
observations. The two-point correlation function of galaxies in the simulation evolves very little since z = 3 and 
has a shape close to a power-law over four orders of magnitude in amplitude. At the present day, the galaxy 
correlation function in the simulation is antibiased relative to the mass on small scales and unbiased on large 
scales. It provides a reasonable match to observations. 

Subject headings: methods: numerical - galaxies: formation: kinematics and dynamics - cosmology: theory - 
hydrodynamical simulation 
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1. INTRODUCTION 

Studies of galaxy formation have advanced at an unprece- 
dented rate in the past few years. Data from the Keck and 
Hubble Space telescopes have revolutionised our view of the 
high redshift Universe (e.g. Steidel et al. 1998) and have led to 
claims that the main phases of galaxy formation activity may 
have now been observed (Baugh et al. 1998). From the the- 
oretical point of view, modelling galaxy formation presents a 
formidable challenge because it involves the synthesis of a wide 
range of disciplines, from early universe cosmology to the mi- 
crophysics and chemistry of star formation. 

Because of the strongly non-linear and asymmetric nature of 
gravitational collapse, the problem of galaxy formation is best 
addressed by direct numerical simulation. The main difficulty 
stems from the huge range of scales spanned by the relevant 
processes, from star formation to large-scale clustering, which 
cannot all be simultaneously resolved with current simulation 
techniques. Two complementary strategies have been devel- 
oped to deal with processes occurring below the resolution limit 
of a simulation. In one of them, a semi-analytic model of the 
dynamics of gas and star formation is used, for example, in con- 
junction with N-body simulations of the formation of dark mat- 
ter halos (Kauffmann et al. 1997, 1999a,b; Benson et al. 1999). 
This technique permits a large dynamic range to be followed, at 
the expense of simplifying assumptions, such as spherical sym- 
metry, for the treatment of the dynamics of cooling gas and star 
formation. The alternative approach, which is the one adopted 
in this paper, is to solve directly the evolution equations for 
gravitationally coupled dark matter and dissipative gas. This 
enables the dynamics of the gas to be treated with a minimum 
of assumptions, at the expense of a severe reduction in the ac- 
cessible dynamic range. As in the semi-analytic approach, a 
phenomenological model for star formation and feedback is re- 
quired. 

Eulerian and Lagrangian numerical hydrodynamics have 
been used to simulate galaxy formation. At present only the 



latter, implemented by means of the Smooth Particle Hydro- 
dynamics or SPH technique, provides sufficient resolution to 
follow the formation of individual galaxies. For example, the 
best Eulerian simulations to date, such as those of Blanton et 
fl/. (1999), have gas resolution elements of ^ 300-500 kpc, 
whereas the early SPH simulation of Carlberg et al. (1990) had 
a spatial resolution of 20 kpc. This, together with the simu- 
lations of Katz et al. (1992), Evrai'd et al. (1994) and Frenk et 
al. (1996), were the first to resolve individual galaxies in rela- 
tively large volumes, allowing detailed studies of the distribu- 
tion of galaxies and the dynamics of galaxies in clusters. How- 
ever, the volumes modelled in this early work were much too 
small to allow reliable investigations of galaxy clustering at the 
present day. Galaxy clustering at high redshift has been inves- 
tigated in the simulations of Katz et al. (1999). 

In this letter we present the first results of a large N- 
body/SPH simulation of galaxy formation in a representative 
volume of a cold dark matter universe, employing about an or- 
der of magnitude more particles than the largest previous study 
of this kind (Katz et al. 1996). Our simulation produced 2266 
galaxies at the present day, compared to 60 in the simulation 
of Katz et al. (1992) and 58 in those of Katz et al. (1996), the 
only other large SPH calculations to have been evolved to the 
present. Earlier dark matter simulations by the Virgo Consor- 
tium (Jenkins et al. 1998) demonstrated the kind of biases re- 
quired for CDM universes to provide a good match to observa- 
tions of galaxy clustering. Here we show that these biases arise 
quite naturally. 

2. THE SIMULATION 

We have simulated a region of a CDM universe with the same 
cosmological parameters as the ACDM simulation of Jenkins 
et al. (1998): mean mass density parameter, VIq = 0.3; cosmo- 
logical constant, A/{3Hq) = 0.7; Hubble constant (in units of 
100 km s"' Mpc"'), h = 0.7 (hereafter we adopt this value, un- 
less otherwise stated), and rms linear fluctuation amplitude in 
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Sh~^ Mpc spheres, erg = 0.9. The baryon fraction was set to 
riijh^ = 0.015, from Big Bang nucleosynthesis constraints (Copi 
et al. 1995). We assumed an unevolving gas metallicity of 0.3 
times the solar value. The simulation was carried out using 
"parallel Hydra", an adaptive, particle -particle, particle-mesh 
N-body/SPH code (Pearce & Couchman 1998), based on the 
publicly released serial version of Couchman et al. (1995). 

Our simulation followed 2097152 dark matter particles and 
the same number of gas particles in a cube of side 100 Mpc and 
required 12492 timesteps to evolve from z = 50 to z = 0. The 
gas mass per particle is 2 x 10^ Mq and, since we typically 
smooth over 32 SPH neighbors, the smallest resolved objects 
have a gas mass of 6.4 x 10'° Mq. We employed a comoving p- 
spline gravitational softening, equivalent to a Plummer soften- 
ing of 14.3 kpc until z=2.5. Thereafter, the softening remained 
fixed at this value, in physical coordinates, and the minimum 
SPH resolution was set to match this value. With our chosen 
parameters, our simulation was able to follow the cooling of 
gas into galactic dark matter halos. The resulting "galaxies" 
typically have 50- 1 000 particles. With a spatial resolution of 
14.3 kpc, we cannot resolve the internal structure of galaxies 
and we must be cautious about the possibility of enhanced tidal 
disruption, drag, and merging within the largest clusters. How- 
ever, as we argue below, there is no evidence that this is a major 
problem. 

As in all studies of this type, a phenomenological model is 
required to treat physical processes occurring below the reso- 
lution limit of the simulation. The first of these is the runaway 
cooling instability present in hierarchical clustering models of 
galaxy formation. At high redshift, the cooling time in dense 
subgalactic objects is so short that most of the gas would cool 
(and presumably turn into stars) unless other processes acted 
to counteract cooling (White & Rees 1978, Cole 1991, White 
& Frenk 1991). Since all the gas in the universe has clearly 
not cooled into dark matter halos, a common assumption is that 
feedback from early generations of stars will have reheated the 
gas, preventing it from cooling catastrophically. 

Although a variety of prescriptions have been used to model 
feedback (e.g. Navarro & White 1993, Steinmetz & Muller 
1995, Katz et al. 1996), this process remains poorly understood. 
In cosmological SPH simulations, gas can only cool efficiently 
in objects above the minimum resolved gas mass, in our case 
6.4 X 1O''^M0. Thus, resolution effects alone act as a crude 
form of feedback. Semi-analytical models of galaxy formation 
suggest that feedback is relatively unimportant on mass scales 
above our resolution limit. We do not, therefore, impose any 
prescription for feedback over and above that provided natu- 
rally by resolution effects. If the rate at which gas cools in the 
simulation is identified with the rate at which stars form, our 
adopted parameters give rise to a cosmic history of star forma- 
tion in broad agreement with data from z ~ 4 to the present 
(Madauef al. 1998). 

The second sub-resolution process that we must model is star 
formation and the associated coupling of different gas phases 
in the interstellar medium. Like feedback, this is a complex 
and poorly understood phenomenon. In some SPH simulations, 
groups of cooled gas particles have been identified with galax- 
ies (the "globs" of Evrard et al. 1994). One disadvantage of this 
procedure is that dense knots of cold gas can affect the cooling 
of surrounding hot material because of the smoothing inherent 
in the SPH technique. To avoid this problem, an alternative 
strategy often used is to assume that gas that has cooled turns 
into colUsionless "stars" according to some heuristic algorithm 
(Navarro & White 1993, Katz et al. 1996, Steinmetz & Muller 
1995). This prescription effectively decouples the cooled gas 
from the hot component. A drawback is that "stellar" systems 



made up of only a few particles are fragile and can easily be 
disrupted. 

We have adopted an intermediate strategy intended as a com- 
promise between the extremes of letting clumps of cool gas 
persist in the simulation and turning them into stars. As in 
the first case, we identify galaxies with groups of gas particles 
that have cooled below lO'^K. However, when computing the 
SPH density of particles with temperatures above lO^K, we do 
not include any contribution from particles below lO^K. All 
other SPH interactions remain unaffected. As in the case where 
cool gas is turned into stars, our model effectively decouples the 
galactic material from the surrounding hot halo gas, but unlike 
this case, "galaxies" in our model are made of dissipative ma- 
terial and thus are more resilient to tidal interactions and merg- 
ers than model stellar galaxies. Our model of the intergalactic 
medium can be regarded as a simple, first step towards a multi- 
phase implementation of SPH, an important requirement when 
dealing with situations in which there are steep density gradi- 
ents. The main effect of our treatment of cool gas is to prevent 
the formation of very massive galaxies in the centers of the rich- 
est clusters in the simulation, as happened, for example, in the 
simulation of Frenk et al. (1996). Thacker et al. (1998) present 
a more detailed discussion of the effects of runaway cooling 
and the production of supermassive objects. 
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Fig, 1. The mass fraction in the tbrm of cold gas in virialized halos, normalized to the 
average baryon fraction in the simulation. The dots show results from the N-body/SPH 
simulation. The pentagons show results from a semi-analytic model constrained to match 
the parameters of the simulation. The open squares show results from a full semi-analytic 
model. 

As a test of our techniques, Fig. 1 compares the amount of 
cold gas in halos in our simulation with that predicted by the 
semi-analytic model of Cole et al. (1999). Halos in the sim- 
ulation were located by first identifying suitable centers using 
the friends-of-friends group finder of Davis et al. (1985) with 
linking parameter, = 0.05, and then growing spheres around 
these centers out to the virial radius (defined as the radius within 
which the mean overdensity is 323; see Eke et al. 1996). Only 
halos with more than 50 dark matter particles were considered, 
of which there are 2353 in the simulation, spanning over 3 or- 
ders of magnitude in mass. The simulation results (shown as 
dots in the figure) are compared with two versions of the semi- 
analytic model. In the first (filled pentagons), the parameters of 
the semi-analytic model were set so as to mimic the conditions 
in the simulation as closely as possible. The mass resolution of 
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the semi-analytic model was degraded to that of the simulation 
and the feedback was switched off. The agreement between 
the simulation and modified semi-analytic model indicates that 
there are only minor differences in the cooling properties of 
the gas calculated with these two different techniques. The 
second comparison is intended as a test of how well the arti- 
ficial "feedback" produced by resolution effects compares with 
the physically motivated feedback prescription used in semi- 
analytic models. In this case (open squares), the parameters of 
the semi-analytic model were set so as to obtain a good match 
to the faint end of the galaxy luminosity function, as discussed 
by Cole et al. (1999). The agreement with the simulation in this 
case is moderate. Clearly, feedback in the semi-analytic model 
prevents cooling within galaxy halos more efficiently than do 
resolution effects in the simulation, but the difference is only 
about 50% for halos similar to that of the Milky Way. 
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Fig. 2. The ratio of gas mass to darlt matter witiiin the viriai radius of each halo, in 
units of the mean baryon fraction in the simulation. The circles show the total gas mass 
fraction while the crosses and triangles show the mass fractions of gas hotter and cooler 
than 12000K respectively. The solid line shows the resolution limit of 32 gas particles. 

3. RESULTS 

The ability of gas to cool is a strong function of the mass 
of the host halo. Fig. 2 shows the fraction of hot and cold gas 
within the viriai radius of each halo, normaUzed to the mean 
baryon fraction of the simulation (10%), as a function of halo 
mass. In small halos just above the resolution limit (indicated 
by the sohd fine), most of the gas cools, but the fraction of cold 
gas decreases rapidly with halo mass as the cooling time for the 
hot gas increases. In large halos, most of the gas never cools. 
The crossover point occurs at a halo mass of ~ 10^^ Mq. Be- 
cause of the generally asymmetric and chaotic nature of halo 
formation, a few low mass halos have baryon fractions in ex- 
cess of the universal mean, but in most galactic halos the baryon 
fraction ranges between 80 and 100% of the cosmic mean. On 
the scale of galaxy clusters, the baryon fraction is 85%, simi- 
lar to the values obtained by White et al. (1993) and Frenk et 
al. (1999). 

We identify "galaxies" in our simulation with dense knots of 
cold gas. These are very easy to locate except in the minority 
of cases where a merger is ongoing or where the galaxy is ex- 
periencing significant tidal disruption or ablation within a clus- 
ter halo. To find galaxies we used the friends-of-friends group 
finder, with a linking length of 0.0164(1 + z) times the mean 



comoving interparticle separation. This selects material with 
overdensity greater than ^ 10^ at z = 0. The galaxy catalogue 
is almost unaffected by large changes in the maximum linking 
length. At the end of the simulation there were 2266 resolved 
objects within the volume. 
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Fig. 3. A comparison between the K-band galaxy luminosity function in the simulation 
with observations. The simulation data are shown by open triangles and the data from 
Gardner et al. (1997) by filled squares. A luminosity normalization factor of T = 2.8 has 
been assumed. Poisson errors are shown. 

We can assign a luminosity to each galaxy in our simula- 
tion using the stellar population synthesis model of Bruzual & 
Chariot (1993). For this purpose we assume that at each model 
output, a fraction, 1/T, of the gas that has cooled since the pre- 
vious output turns into stars in an instantaneous burst with a 
Salpeter IMF whose subsequent spectrophotometric evolution 
is given by the synthesis model. Because the output times were 
relatively infrequent, this procedure works best for K-band lu- 
minosities. 

In Fig. 3 we compare the resulting K-band galaxy luminosity 
function with the observational data of Gardner et al. (1997). 
The shape of the model luminosity function agrees well with 
the data, and the model and observed functions match well if 
we set the luminosity normahzation factor, T = 2.8. This im- 
plies that only 35% of the cold gas has been turned into visible 
stars, with the rest remaining in dense gas clouds and brown 
dwarfs or hidden in some other form. The associated mass-to- 
light ratios are about twice as large as those measured for el- 
liptical galaxies, but these numbers should be treated with cau- 
tion. Not only is our star formation prescription very crude, 
but our model ignores the effects of metallicity and obscuration 
by dust. Furthermore, as the comparison with the full semi- 
analytic model indicates, too much gas has probably cooled in 
our simulation because of our neglect of feedback processes. In 
spite of these reservations, the agreement in Fig. 3 is very good 
and suggests that our simulation provides a realistic description 
of the formation of bright galaxies. 

The relatively large volume of our simulation allows a reU- 
able measurement of the clustering properties of galaxies and 
their relation to the clustering properties of the mass. The 
galaxy and mass two-point correlation functions at various 
epochs are plotted in Fig. 4. The mass correlation function 
agrees very well with the results of our earlier, dark matter only 
simulations which followed a cubic region of side 342 Mpc us- 
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ing 16.8 million particles (Jenkins et al. 1998). The clustering 
amplitude of the mass grows by a factor of about 30 between 
the two epochs shown in the figure, z = 3 and z = 0. By contrast, 
the galaxy correlation function hardly evolves at all between 
z = 3 and z = 0. 
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Fig. 4. Mass and galaxy correlation functions. The dashed lines show the mass cor- 
relation functions at z = 3 and z = 0. The soHd, long dashed and dotted Unes show the 
galaxy correlation functions in the simulation at the indicated redshifts. The squares show 
the observed, real-space correlation function, estimated by Baugh (1996), from the APM 
survey. 

The difference between the clustering growth rates of galax- 
ies and mass is a manifestation of "biased galaxy formation", 
the preferential formation of galaxies in high peaks of the pri- 
mordial density field. It was already seen in the first simulations 
of cold dark matter models by Davis et al. (1985), in which 
galaxies were put in "by hand" near high peaks of the initial 
density field. It is also very clear in the SPH simulations of 
Evrard et al. (1994) and Katz et al. (1999), and can even be in- 
ferred from the N-body only simulations of Bagla (1998) and 
Colin et al. (1999). Semi-analytic techniques have been used 
on their own (Baugh et al. 1999), or combined with N-body 
simulations (Kauffmann et al. 1999b), for detailed study of the 
clustering evolution of galaxies, while specific applications to 
high redshift Lyman-break galaxies are to be found in Baugh et 



al. (1998) and Governato et al. (1998). The latter models pro- 
vide an excellent description of the strong clustering discovered 
by Adelberger et al. (1998). 

In Fig. 4 we also plot the observed, real-space galaxy corre- 
lation function at z — 0, estimated by Baugh (1996) from the 
APM survey (filled squares.) This may be compared with the 
z = results in our simulation (solid line). On scales larger than 
a few hundred kpc the agreement is good. (The differences at 
r > I0h~^ Mpc are due, for the most part, to finite volume ef- 
fects, as we have verified by comparison with the larger sim- 
ulations of Jenkins et al. 1998.) Beyond l/z~' Mpc the galaxy 
correlation function is very close to the mass correlation func- 
tion. On smaller scales galaxies are less strongly clustered than 
the mass, or antibiased, an effect that persists until separations 
of ^ 100/!~'kpc. At small separations the model correlations lie 
above the APM data. Over nearly four orders of magnitude in 
amplitude, the model galaxy correlation function is very close 
to a power-law even though the mass correlation function is not. 
An essentially featureless galaxy correlation function was also 
obtained for the same cosmological model in the semi-analytic 
model of Benson et al. (1999) and, for some parameter combi- 
nations, in those of Kauffmann et al. (1999a). 

4. CONCLUSIONS 

The simulation presented here is the first to resolve galaxy 
formation in a large enough volume to allow a reliable study 
of the demographics and clustering of galaxies. Our results are 
encouraging: the resulting luminosity function and correlation 
function of galaxies are broadly consistent with observations. 
Furthermore, the correlation function of bright galaxies in the 
simulation changes little since z = 3, in agreement with results 
from semi-analytic studies and with available data at high red- 
shift. Further progress will require a more detailed treatment 
of the astrophysics of galaxy formation, particularly of the pro- 
cesses of star formation and feedback. 
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